Examining site effects in DTIMSCAMRAS
Virgilio Gonzenbach
2023-06-06
library(tidyverse)
library(lme4)
options(scipen=999)
df = read.csv(here::here('avg_tensor_by_roi_wide.csv'),
colClasses = c('subject' = 'factor',
'site' = 'factor',
'visit' = 'factor'))
outcomes = df %>%
select(where(is.numeric)) %>%
colnamesSummary
The purpose of the DTIMSCAMRAS study is to test for site effects on DTI metrics of different brain regions in the context of a traveling subjects design with a harmonized imaging protocol. Participants (n=11) with stable Multiple Sclerosis (MS) were scanned in 4 different locations: National Institute of Health (NIH), Brigham and Women’s Hospital (BWH), Johns Hopkins University (Hopkins), and The University of Pennsylvania (Penn). Two scans were collected per site visit, and participants were re-positioned between scans. Penn, BWH and NIH used Siemens scanners, while Hopkins used a Phillips scanner. The imaging modalities that were collected include T1s (MPRAGE), FLAIR, and Diffusion-weighted images (DWI).
Preprocessing
To preprocess imaging data, qsiprep
was used for bias correction, skull-stripping and co-registration of the
T1s and DWIs, and specialized denoising and artifact correction for the
DWIs. Then, DTIFIT was used to compute the following scalar
maps:
- Fractional Anisotropy
- Mean Diffusivity
- Axial Diffusivity
- Radial Diffusivity
Concurrently, several segmentation pipelines were used on the T1s (+ FLAIRs in the case of MIMOSA) with the purpose of defining particular regions of interest (ROIs):
- White Matter and Gray Matter Segmentation methods:
- ANTs ATROPOS
- FSL FAST
- Multi-Atlas Segmentation with Joint Label Fusion using the MUSE templates (JLF WM and JLF GM)
- Thalamus Segmentation
- FSL FIRST
- Multi-Atlas Segmentation with Joint Label Fusion using the OASIS templates (JLF THALAMUS)
- Lesion Segmentation methods:
- MIMOSA
To remove lesions appearing in GM or WM (or thalamus) from the calculation, mimosa masks were inverted and applied to the segmentation label maps in order to zero out the tissue masks at the location of each lesion. These adjusted segmentation results were then used to average the DTI scalar maps across all voxels belonging to each of the ROIs. The 36 resulting DTI metrics make up the set of outcome variables which are the focus of the site effects test. Below, the outcome variables are listed.
data.frame(OUTCOME = outcomes) %>%
separate(OUTCOME, c("SEGMENTATION", "ROI", "SCALAR"), remove = FALSE)Note the general
Testing site effects
Testing for site effects was performed using two approaches: a permutation-based approach and a deviance test of mixed models.
Additionally, as Hopkins data was acquired using a Philips scanner, which is in contrast to all other sites, the two testing approaches were performed again while excluding Hopkins data.
Permutation-Based F-test
A permutation-based F-test was used to test for site effects. For each subject, the squared difference between all possible pairs of intra-site and inter-site measures were computed. For instance, take the first row of the subject data.frame below, which has BWH as the “reference site”: the squared differences between that row and the 6 non-BWH rows are computed for \(y\)——this process is repeated for all rows (and the results averaged) to arrive to the mean inter-site squared difference for this subject; to calculate the mean intra-site squared differences, the squared differences of the 4 intrasite pairs are averaged.
df %>%
filter(subject == '01001') %>%
select(!where(is.numeric), ATROPOS_GM_AD) %>%
rename(y = ATROPOS_GM_AD)These differences were averaged across all subjects resulting in the Mean Squared difference Between Sites (\(\overline{MS_{B}}\)) and the Mean Squared difference Within Sites (\(\overline{MS_{W}}\)). The test statistic is composed of their ratio:
\[ F = \frac{\overline{MS_{B}}}{\overline{MS_{W}}} \]
For each metric, the observed test statistic was compared to a distribution of test statistics derived from 10000 permutations (i.e. a null distribution). The site labels were permuted within each subject. The proportion of null results higher than the observed test statistic served as a preliminary p-value \(p_0\). To prevent p-values of 0 (when the observed statistic was higher than all permuted results), the p-value was shifted using the following formula:
\[ p = \frac{10000*p_0 + 1}{10000 + 1}\]
Crossed-design Mixed Linear Models
As an additional test of site effects, mixed models were performed
corresponding to a crossed design using the lmer function
of the lme4 package. For each outcome variable, two
mixed-effects models were fitted: The formula for the base model
includes a random intercept for site as well as a random intercept for
subject. The formula for the extended model also includes these terms in
addition to a fixed effect term for site.
Base Model:
y ~ (1|subject) + (1|site)Extended model:
y ~ site + (1|subject) + (1|site)For each outcome, these two models were compared through a deviance
test using the anova function in R. P-values for these
tests are reported as a test of site effects.
Source code for this project can be found on GitHub.
Key Findings
- Permutation tests revealed significant site effects in 31 out of 36 outcomes after FDR adjustment when all data is included. After excluding Hopkins data, the number of significant site effects decreased to 13. Most of the DTI metrics (i.e. outcomes) that showed differences were derived from thalamus segmentation results.
- For mixed models, most outcome variables (13 out of 36) show a significant test after FDR adjustment indicating the presence of site effects when all data is included. When Hopkins data is excluded, none of the outcome variables show site effects.
Questions for follow up
- It is unclear to what extent potential segmentation differences across sites might contribute to site effects, i.e. to what extent do site differences stem from ROI coverage differences (due to scanner effects on the T1s), and to what extent are the differences due to scanner effects on DWI images?
Data Inventory
Subjects 03001 and 03002 did not complete imaging at all 4 sites. Subject 03001 completed imaging at Penn and BWH only and 03002 completed imaging at BWH, Hopkins, and the NIH only.
Subject 02001 is missing DTI data from their NIH visit.
Scans at site per subject
t_df <- df %>%
count(subject, site, .drop = FALSE) %>%
tidyr::pivot_wider(names_from = site, values_from = n) %>%
mutate(`Row Totals` = BWH + Hopkins + NIH + Penn)
t_df <- t_df %>%
summarize(across(where(is.integer), sum)) %>%
bind_rows(t_df, .) %>%
mutate(across(where(is.factor), as.character))
t_df[is.na(t_df)] <- "Column Totals"
t_dfFor Atropos (substracting failed segmentations)
Atropos segmentation failed for 6 images: 04001NIH01, 01003NIH01, 03002NIH02, 04003NIH01, 04001BWH02, and 03001BWH01.
knitr::include_graphics(here::here(c('misc/atropos_success.png', 'misc/atropos_fail.png')))
Left: Example of successful Atropos segmentation; Right: Example of a failed Atropos segmentation with one ROI missing
The following table excludes the failed atropos segmentations. Note subject 03001 only has 3 images after excluding failed atropos segmentations.
t_df <- df %>%
unite(sub, subject, site, visit, sep = '-') %>%
filter(!sub %in% c('04001-NIH-01', '01003-NIH-01', '03002-NIH-02', '04003-NIH-01', '04001-BWH-02', '03001-BWH-01')) %>%
separate(sub, c('subject', 'site', 'visit')) %>%
mutate_if(is.character, as.factor) %>%
count(subject, site, .drop = FALSE) %>%
tidyr::pivot_wider(names_from = site, values_from = n) %>%
mutate(`Row Totals` = BWH + Hopkins + NIH + Penn)
t_df <- t_df %>%
summarize(across(where(is.integer), sum)) %>%
bind_rows(t_df, .) %>%
mutate(across(where(is.factor), as.character))
t_df[is.na(t_df)] <- "Column Totals"
t_dfThe failed atropos segmentations are included in the following visualizations but were excluded for the site effects tests.
TODO: compare harmonization with and without hopkins
Visualizations
Raw Data
Fractional Anisotropy
Densities
Densities are colored by site, while the black line represents the overall density aggregated across sites. Distribution of the different sites tend to cluster together except for JLF Thalamus. Additionally, note the long tails caused by outliers in the ATROPOS distributions.
plot_densities <- function(var){
df %>%
ggplot(aes_string(x=var, color='site')) +
geom_density() +
stat_density(aes_string(x = var), geom = "line", color = "black") +
theme_bw() +
xlab(str_replace_all(var, "_", " "))
}
vars <-df %>%
select(ends_with('FA')) %>%
colnames()
dplots <- lapply(vars, plot_densities) %>%
setNames(vars)for (name in names(dplots)){
cat("#####", str_replace_all(name, "_", " "), "\n")
print(dplots[[name]])
cat("\n\n")
}ATROPOS GM FA
ATROPOS WM FA
FAST GM FA
FAST WM FA
FIRST THALAMUS FA
JLF GM FA
JLF WM FA
JLF THALAMUS FA
MIMOSA LESION FA
Box plots
Box plots show difference in medians across sites for most
variables.
Additionally, note the severe outliers in the ATROPOS metrics which
correspond to the failed segmentations. For FIRST THALAMUS, subject
01002 could be considered an outlier; their FA values are consistent
across sites and scans, however, suggesting this variation is meaningful
and that this subject should be retained.
plot_avg_tensor_values <- function(tensor){
df %>%
select(!is.numeric, ends_with(tensor)) %>%
pivot_longer(ends_with(tensor), names_to = 'seg', values_to = 'average') %>%
separate(seg, into = c('segmentation', 'roi', 'tensor')) %>%
unite(segmentation, segmentation, roi, sep = " ") %>%
mutate(segmentation = factor(segmentation,
levels = c("ATROPOS GM", "ATROPOS WM", "FAST GM", "FAST WM",
"JLF GM", "JLF WM", "JLF THALAMUS", "FIRST THALAMUS", "MIMOSA LESION"))) %>%
ggplot(aes(x=site, y=average)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(aes(color = subject, shape=visit), alpha=0.5) +
#facet_grid(site~.) +
facet_grid(segmentation ~ .) +
coord_flip() +
ggtitle(sprintf("Mean %s in ROI by site", tensor)) +
theme_bw() +
# theme(axis.text.x = element_text(angle = 60, vjust = 1, hjust=1)) +
xlab("") +
ylab("") +
theme(strip.text.y = element_text(angle = 0)) +
scale_x_discrete(expand = c(0.15, 0.15))
}
plot_avg_tensor_values('FA')
Mean Diffusivity
Densities
For Mean Diffusivity, distributions vary more widely across sites. In particular, MIMOSA LESION MD values are quite different for Hopkins data.
vars <-df %>%
select(ends_with('MD')) %>%
colnames()
dplots <- lapply(vars, plot_densities) %>%
setNames(vars)for (name in names(dplots)){
cat("#####", str_replace_all(name, "_", " "), "\n")
print(dplots[[name]])
cat("\n\n")
}ATROPOS GM MD
ATROPOS WM MD
FAST GM MD
FAST WM MD
FIRST THALAMUS MD
JLF GM MD
JLF WM MD
JLF THALAMUS MD
MIMOSA LESION MD
Boxplot
As before, medians appear to be different across sites.
plot_avg_tensor_values('MD')
Radial Diffusivity
Densities
Radial Diffusivity also shows some variations across sites (perhaps less than MD). JLF GM RD in NIH shows multiple peaks, and JLF WM RD shows skewed distributions for all sites. As before the metric for MIMOSA shows a different distribution for Hopkins data.
vars <-df %>%
select(ends_with('RD')) %>%
colnames()
dplots <- lapply(vars, plot_densities) %>%
setNames(vars)for (name in names(dplots)){
cat("#####", str_replace_all(name, "_", " "), "\n")
print(dplots[[name]])
cat("\n\n")
}ATROPOS GM RD
ATROPOS WM RD
FAST GM RD
FAST WM RD
FIRST THALAMUS RD
JLF GM RD
JLF WM RD
JLF THALAMUS RD
MIMOSA LESION RD
Boxplot
As before, medians appear to be differet across sites.
plot_avg_tensor_values('RD')